Spatial Ecology and Macroecology

Practical - Week 2

Florencia Grattarola, Friederike Wölke, Gabriel Ortega

(Department of Spatial Sciences)

2023-10-09

What are we going to see today?

How to plot biodiversity data, explore patterns at different resolutions, and make pretty maps.

  1. Create a grid over the study area (extent/grain)
  2. Calculate biodiversity metrics per area (species richness, endemism)
  3. Visualisation/mapping (color schemes)
  4. Description of patterns (hotspots)

Exercises:

  1. LOCAL scale - occurrence records 🦌
  2. GLOBAL scale - expert range maps 🐢

Some preparation before starting to code

Data download

Please download the following files and store them on you data folder:

  • mammalsCZ.rds
  • CZE_adm0.gpkg
  • KvadratyCR_JTSK.gpkg
  • testudines.gpkg

Spatial analyses

For geocomputation today we will use the sf package.

install.packages('sf') # install the package if you haven't already
library(sf) # load
Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
packageVersion('sf')
[1] '1.0.8'


sf_use_s2(FALSE) # switch spherical geometry off
Spherical geometry (s2) switched off

Spatial data download

For downloading country-specific and world polygons we will use the rnaturalearth package.

install.packages(c('rnaturalearth', 'rnaturalearthdata')) # the second package is also needed
library(rnaturalearth) # load
packageVersion('rnaturalearth')
[1] '0.3.2'

Visualisation

For all the map plots today we will use the ggplot package. It’s part of the tidyverse.

library(tidyverse) # load

LOCAL Mammal’s of the Czech Republic

Photo of Capreolus capreolus observed in Czech Republic by romanvrbicek licensed CC-BY-NC, via iNaturalist

Mapping occurrence records 🦌

Steps:

  1. Get occurrence records of mammals in the Czech Republic (POINTS)
  2. Create a grid layer of the Czech Republic (POLYGONS)
  3. Calculate sampling effort (N) and species richness per grid cell (SR)
  4. Visualise N and SR patterns in Czech Republic

Mapping occurrence records 🦌

  1. Get occurrence records of mammals in the Czech Republic

Get the data from last practical (downloaded from GBIF)

mammalsCZ <- readRDS('data/mammalsCZ.rds')
mammalsCZ
# A tibble: 3,931 × 177
   key    scien…¹ decim…² decim…³ issues datas…⁴ publi…⁵ insta…⁶ hosti…⁷ publi…⁸
   <chr>  <chr>     <dbl>   <dbl> <chr>  <chr>   <chr>   <chr>   <chr>   <chr>  
 1 40115… Dama d…    49.2    16.5 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 2 40116… Castor…    50.2    14.6 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 3 40150… Myocas…    49.7    15.1 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 4 40181… Myocas…    50.1    14.4 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 5 40149… Sus sc…    49.2    16.5 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 6 40149… Dama d…    49.2    16.5 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 7 40149… Capreo…    49.6    16.7 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 8 40149… Lepus …    49.6    16.7 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
 9 40149… Myocas…    50.1    14.6 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
10 40148… Myocas…    49.8    14.7 cdc,c… 50c950… 28eb1a… 997448… 28eb1a… US     
# … with 3,921 more rows, 167 more variables: protocol <chr>,
#   lastCrawled <chr>, lastParsed <chr>, crawlId <int>, basisOfRecord <chr>,
#   occurrenceStatus <chr>, taxonKey <int>, kingdomKey <int>, phylumKey <int>,
#   classKey <int>, orderKey <int>, familyKey <int>, genusKey <int>,
#   speciesKey <int>, acceptedTaxonKey <int>, acceptedScientificName <chr>,
#   kingdom <chr>, phylum <chr>, order <chr>, family <chr>, genus <chr>,
#   species <chr>, genericName <chr>, specificEpithet <chr>, taxonRank <chr>, …

Mapping occurrence records 🦌

  1. Get occurrence records of mammals in the Czech Republic

Let’s keep only a few fields

mammalsCZ <- mammalsCZ %>% select(species, order, eventDate, 
                                  decimalLongitude, decimalLatitude,
                                  stateProvince, countryCode, datasetName)

Mapping occurrence records 🦌

  1. Get occurrence records of mammals in the Czech Republic

We will transform our data (table) into an sf object using st_as_sf()

st_as_sf(
  x,
  ...,
  agr = NA_agr_,
  coords,
  wkt,
  dim = "XYZ",
  remove = TRUE,
  na.fail = TRUE,
  sf_column_name = NULL
)

Mapping occurrence records 🦌

  1. Get occurrence records of mammals in the Czech Republic
mammalsCZ %>% 
  filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>% # filter records without coordinates
  st_as_sf(coords=c('decimalLongitude', 'decimalLatitude'),
           crs=4326)
Simple feature collection with 3930 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 12.20682 ymin: 48.63992 xmax: 18.79605 ymax: 50.99892
Geodetic CRS:  WGS 84
# A tibble: 3,930 × 7
   species       order event…¹ state…² count…³ datas…⁴            geometry
 * <chr>         <chr> <chr>   <chr>   <chr>   <chr>           <POINT [°]>
 1 Dama dama     Arti… 2023-0… Jihomo… CZ      iNatur… (16.52097 49.19989)
 2 Castor fiber  Rode… 2023-0… Středo… CZ      iNatur… (14.64081 50.21619)
 3 Myocastor co… Rode… 2023-0… Středo… CZ      iNatur… (15.08824 49.73967)
 4 Myocastor co… Rode… 2023-0… Prague  CZ      iNatur… (14.41362 50.08528)
 5 Sus scrofa    Arti… 2023-0… Jihomo… CZ      iNatur… (16.54771 49.19072)
 6 Dama dama     Arti… 2023-0… Jihomo… CZ      iNatur… (16.54895 49.19099)
 7 Capreolus ca… Arti… 2023-0… Pardub… CZ      iNatur… (16.72316 49.63749)
 8 Lepus europa… Lago… 2023-0… Pardub… CZ      iNatur… (16.72331 49.63835)
 9 Myocastor co… Rode… 2023-0… Prague  CZ      iNatur… (14.60999 50.07419)
10 Myocastor co… Rode… 2023-0… Středo… CZ      iNatur… (14.65434 49.78069)
# … with 3,920 more rows, and abbreviated variable names ¹​eventDate,
#   ²​stateProvince, ³​countryCode, ⁴​datasetName

Note that crs=4326 is related to the EPSG:4326 CRS: WGS 84

Mapping occurrence records 🦌

  1. Get occurrence records of mammals in the Czech Republic

And we save

mammalsCZ_sf <- mammalsCZ %>% 
  filter(!is.na(decimalLongitude) & !is.na(decimalLatitude)) %>% # filter records without coordinates
  st_as_sf(coords=c('decimalLongitude', 'decimalLatitude'),
           crs=4326) # this is the projection that the data is in

Mapping occurrence records 🦌

  1. Create a grid layer of the Czech Republic

Load the administrative borders using st_read

CZ_borders <- st_read('data/CZE_adm0.gpkg')
Reading layer `CZE_adm0' from data source 
  `/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/CZE_adm0.gpkg' 
  using driver `GPKG'
Simple feature collection with 1 feature and 70 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 12.08586 ymin: 48.54084 xmax: 18.86253 ymax: 51.05438
Geodetic CRS:  WGS 84

We’ll see how to get this polygon data directly from R later

Mapping occurrence records 🦌

  1. Create a grid layer of the Czech Republic
CZ_borders
Simple feature collection with 1 feature and 70 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 12.08586 ymin: 48.54084 xmax: 18.86253 ymax: 51.05438
Geodetic CRS:  WGS 84
  GADMID ISO     NAME_ENGLI       NAME_ISO       NAME_FAO      NAME_LOCAL
1     59 CZE Czech Republic CZECH REPUBLIC Czech Republic Ceská republika
       NAME_OBSOL NAME_VARIA NAME_NONLA         NAME_FRENC      NAME_SPANI
1 Moravia|Bohemia      Cesko       <NA> République Tchèque República Checa
  NAME_RUSSI         NAME_ARABI NAME_CHINE      WASPARTOF CONTAINS
1      ????? ????????? ????????      ????? Czechoslovakia     <NA>
       SOVEREIGN ISO2  WWW FIPS ISON  VALIDFR VALIDTO AndyID  POP2000     SQKM
1 Czech Republic   CZ <NA>   EZ  203 19930101 Present     65 10271830 78495.16
   POPSQKM      UNREGION1 UNREGION2 DEVELOPING CIS Transition OECD
1 130.8594 Eastern Europe    Europe          2   0          0    1
               WBREGION            WBINCOME        WBDEBT WBOTHER CEEAC CEMAC
1 Europe & Central Asia Upper middle income Less indebted    <NA>     0     0
  CEPLG COMESA EAC ECOWAS IGAD IOC MRU SACU UEMOA UMA PALOP PARTA CACM EurAsEC
1     0      0   0      0    0   0   0    0     0   0     0     0    0       0
  Agadir SAARC ASEAN NAFTA GCC CSN CARICOM EU CAN ACP Landlocked AOSIS SIDS
1      0     0     0     0   0   0       0  1   0   0          1     0    0
  Islands LDC Shape_Leng Shape_Area                           geom
1       0   0   26.21557   9.813959 MULTIPOLYGON (((13.17648 49...

Mapping occurrence records 🦌

  1. Create a grid layer of the Czech Republic

We will actually load a Czech standard grid layer of 100 km2 area

CZ_grids <- st_read('data/KvadratyCR_JTSK.gpkg')
Reading layer `KvadratyCR_JTSK' from data source 
  `/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/KvadratyCR_JTSK.gpkg' 
  using driver `GPKG'
Simple feature collection with 678 features and 10 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North

We’ll see how to create this grid layer directly from R later

Mapping occurrence records 🦌

  1. Create a grid layer of the Czech Republic
CZ_grids
Simple feature collection with 678 features and 10 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
First 10 features:
   OBJECTID ENTITY    AREA PERIMETE KVADRAT         X         Y Shape_Leng
1         1      3 130.032   45.626    4951 -740241.8 -935317.5   45626.41
2         2      5 130.036   45.627    4952 -728665.0 -936923.5   45627.17
3         3      6 130.303   45.675    5051 -741782.9 -946335.6   45675.21
4         4      8 130.040   45.628    4953 -717084.4 -938503.6   45627.96
5         5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
6         6     12 130.311   45.677    5053 -718576.3 -949528.8   45676.59
7         7     14 130.050   45.630    4955 -693912.1 -941586.4   45629.61
8         8     17 130.055   45.630    4956 -682320.5 -943089.1   45630.48
9         9     18 130.319   45.678    5055 -695354.9 -952618.5   45678.09
10       10     19 130.060   45.631    4957 -670725.5 -944566.0   45631.36
   Shape_Area   ID                           geom
1   130031546 4951 POLYGON ((-745252.4 -928997...
2   130035855 4952 POLYGON ((-733689.7 -930614...
3   130302745 5051 POLYGON ((-746805.7 -940014...
4   130040349 4953 POLYGON ((-722123.3 -932206...
5   130306551 5052 POLYGON ((-735218.5 -941635...
6   130310574 5053 POLYGON ((-723627.4 -943229...
7   130049749 4955 POLYGON ((-698979.2 -935311...
8   130054715 4956 POLYGON ((-687401.8 -936825...
9   130319128 5055 POLYGON ((-700434.2 -946342...
10  130059749 4957 POLYGON ((-675820.7 -938313...

Mapping occurrence records 🦌

  1. Create a grid layer of the Czech Republic

Check the layer’s Coordinate Reference System (CRS)

st_crs(CZ_borders) == st_crs(CZ_grids)
[1] FALSE

Mapping occurrence records 🦌

  1. Create a grid layer of the Czech Republic

Check the layer’s Coordinate Reference System (CRS)

st_crs(CZ_borders)
Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
st_crs(CZ_grids)
Coordinate Reference System:
  User input: S-JTSK / Krovak East North 
  wkt:
PROJCRS["S-JTSK / Krovak East North",
    BASEGEOGCRS["S-JTSK",
        DATUM["System of the Unified Trigonometrical Cadastral Network",
            ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4156]],
    CONVERSION["Krovak East North (Greenwich)",
        METHOD["Krovak (North Orientated)",
            ID["EPSG",1041]],
        PARAMETER["Latitude of projection centre",49.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8811]],
        PARAMETER["Longitude of origin",24.8333333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8833]],
        PARAMETER["Co-latitude of cone axis",30.2881397527778,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",1036]],
        PARAMETER["Latitude of pseudo standard parallel",78.5,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8818]],
        PARAMETER["Scale factor on pseudo standard parallel",0.9999,
            SCALEUNIT["unity",1],
            ID["EPSG",8819]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["GIS."],
        AREA["Czechia; Slovakia."],
        BBOX[47.73,12.09,51.06,22.56]],
    ID["EPSG",5514]]

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)

Frist, we will transform the layer’s CRS to S-JTSK / Krovak East North using st_transform().

mammalsCZ_sf <- st_transform(mammalsCZ_sf,  
                             crs = st_crs(CZ_grids)) # the same CRS as the CZ_grids layer

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)
mammalsCZ_sf
Simple feature collection with 3930 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -896937 ymin: -1224642 xmax: -435147.4 ymax: -943016.8
Projected CRS: S-JTSK / Krovak East North
# A tibble: 3,930 × 7
   species       order event…¹ state…² count…³ datas…⁴             geometry
 * <chr>         <chr> <chr>   <chr>   <chr>   <chr>            <POINT [m]>
 1 Dama dama     Arti… 2023-0… Jihomo… CZ      iNatur… (-604409.5 -1159540)
 2 Castor fiber  Rode… 2023-0… Středo… CZ      iNatur…   (-725332 -1030943)
 3 Myocastor co… Rode… 2023-0… Středo… CZ      iNatur… (-700429.6 -1087674)
 4 Myocastor co… Rode… 2023-0… Prague  CZ      iNatur… (-743384.1 -1043179)
 5 Sus scrofa    Arti… 2023-0… Jihomo… CZ      iNatur… (-602582.8 -1160764)
 6 Dama dama     Arti… 2023-0… Jihomo… CZ      iNatur… (-602489.9 -1160744)
 7 Capreolus ca… Arti… 2023-0… Pardub… CZ      iNatur… (-584603.6 -1112730)
 8 Lepus europa… Lago… 2023-0… Pardub… CZ      iNatur… (-584582.6 -1112636)
 9 Myocastor co… Rode… 2023-0… Prague  CZ      iNatur…   (-729625 -1046301)
10 Myocastor co… Rode… 2023-0… Středo… CZ      iNatur… (-730828.3 -1079075)
# … with 3,920 more rows, and abbreviated variable names ¹​eventDate,
#   ²​stateProvince, ³​countryCode, ⁴​datasetName

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)

Let’s plot the sf objects

Mapping occurrence records 🦌

ggplot() + 
  geom_sf(data=CZ_borders, fill='white') +
  geom_sf(data=CZ_grids, fill=NA) + 
  geom_sf(data=mammalsCZ_sf)

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)

To calculate number of records and number of species per grid cell, we will use st_join()

st_join(
  x,    #   object of class sf
  y,    #   object of class sf
  join = st_intersects,
  ...,
  suffix = c(".x", ".y"),
  left = TRUE,
  largest = FALSE
)

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)

To calculate number of records and number of species per grid cell, we will use st_join()

st_join(CZ_grids,       # POLYGONS (grid cells)
        mammalsCZ_sf)   # POINTS
Simple feature collection with 4183 features and 16 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
First 10 features:
    OBJECTID ENTITY    AREA PERIMETE KVADRAT         X         Y Shape_Leng
1          1      3 130.032   45.626    4951 -740241.8 -935317.5   45626.41
2          2      5 130.036   45.627    4952 -728665.0 -936923.5   45627.17
3          3      6 130.303   45.675    5051 -741782.9 -946335.6   45675.21
4          4      8 130.040   45.628    4953 -717084.4 -938503.6   45627.96
5          5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.1        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.2        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.3        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.4        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
5.5        5      9 130.307   45.676    5052 -730181.5 -947945.1   45675.88
    Shape_Area   ID             species        order           eventDate
1    130031546 4951                <NA>         <NA>                <NA>
2    130035855 4952                <NA>         <NA>                <NA>
3    130302745 5051                <NA>         <NA>                <NA>
4    130040349 4953                <NA>         <NA>                <NA>
5    130306551 5052 Capreolus capreolus Artiodactyla 2020-06-17T00:00:00
5.1  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-20T00:00:00
5.2  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-20T00:00:00
5.3  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-18T00:00:00
5.4  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-23T00:00:00
5.5  130306551 5052 Capreolus capreolus Artiodactyla 2020-06-23T00:00:00
    stateProvince countryCode datasetName                           geom
1            <NA>        <NA>        <NA> POLYGON ((-745252.4 -928997...
2            <NA>        <NA>        <NA> POLYGON ((-733689.7 -930614...
3            <NA>        <NA>        <NA> POLYGON ((-746805.7 -940014...
4            <NA>        <NA>        <NA> POLYGON ((-722123.3 -932206...
5            <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.1          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.2          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.3          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.4          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...
5.5          <NA>          CZ        <NA> POLYGON ((-735218.5 -941635...

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)

Now, we need to summarise the data per grid-cell.
Luckily, tidyverse methods also work for sf objects :)


We will do it with group_by() and summarise().

First, we will group by grid-cells and then we count values per grid.

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)

Let’s calculate the number of records (N) and the number of species per grid-cell (SR)

st_join(CZ_grids, mammalsCZ_sf) %>% 
  group_by(OBJECTID) %>% # the name of the column that has the index
  summarise(N=sum(!is.na(species)), # calculates the number of points in the polygon
            SR=n_distinct(species, na.rm = TRUE))  # calculates the number of different 'species' in the polygon
Simple feature collection with 678 features and 3 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -912065.7 ymin: -1229672 xmax: -421278.2 ymax: -928997.8
Projected CRS: S-JTSK / Krovak East North
# A tibble: 678 × 4
   OBJECTID     N    SR                                                     geom
      <int> <int> <int>                                            <POLYGON [m]>
 1        1     0     0 ((-745252.4 -928997.8, -733689.7 -930614.9, -733689.7 -…
 2        2     0     0 ((-733689.7 -930614.9, -722123.3 -932206.2, -722123.3 -…
 3        3     0     0 ((-746805.7 -940014.3, -735218.5 -941635, -735218.5 -94…
 4        4     0     0 ((-722123.3 -932206.2, -710553.1 -933771.6, -710553.1 -…
 5        5    27     3 ((-735218.5 -941635, -723627.4 -943229.9, -723627.4 -94…
 6        6     5     3 ((-723627.4 -943229.9, -712032.7 -944798.9, -712032.7 -…
 7        7     0     0 ((-698979.2 -935311.3, -687401.8 -936825.2, -687401.8 -…
 8        8     0     0 ((-687401.8 -936825.2, -675820.7 -938313.3, -675820.7 -…
 9        9     0     0 ((-700434.2 -946342, -688832.2 -947859.3, -688832.2 -94…
10       10     0     0 ((-675820.7 -938313.3, -664236.1 -939775.7, -664236.1 -…
# … with 668 more rows

Mapping occurrence records 🦌

  1. Calculate sampling effort (N) and species richness per grid cell (SR)

Let’s store the output into a new object CZ_mammals_grids.

mammalsCZ_grids <- st_join(CZ_grids, mammalsCZ_sf) %>% 
  group_by(OBJECTID) %>% 
  summarise(N=sum(!is.na(species)), 
            SR=n_distinct(species, na.rm = TRUE))

Mapping occurrence records 🦌

  1. Visualise N and SR patterns in Czech Republic

Finally, let’s plot this.
We will do it using geom_sf(), a ggplot function to visualise sf objects

Mapping occurrence records 🦌

  1. Visualise N and SR patterns in Czech Republic
ggplot() + 
  geom_sf(data=mammalsCZ_grids) +
  coord_sf(crs=4326)

Mapping occurrence records 🦌

  1. Visualise N and SR patterns in Czech Republic

Where are the nice colours?

We need to indicate which column from the object should the grids be filled with.

Mapping occurrence records 🦌

  1. Visualise N and SR patterns in Czech Republic
ggplot() + 
  geom_sf(data=mammalsCZ_grids, aes(fill=SR)) +
  coord_sf(crs=4326)

Mapping occurrence records 🦌

  1. Visualise N and SR patterns in Czech Republic

Let’s get a nicer color scale.

Bear in mind that we see colors differently. Thus, it’s important to consider colorblind safe palettes

Check out The R Graph Gallery for many cool graphic options with ggplot

Mapping occurrence records 🦌

  1. Visualise N and SR patterns in Czech Republic
ggplot() + 
  geom_sf(data=mammalsCZ_grids, aes(fill=SR)) +
  scale_fill_fermenter(palette ='YlOrBr', n.breaks=9, direction = 1) + # fill of the grids
  geom_sf(data=CZ_borders, fill=NA) +
  coord_sf(crs=4326) +
  theme_bw()

Now, here’s a better plot :)

Mapping occurrence records 🦌

  1. Visualise N and SR patterns in Czech Republic

The hotspots of species richness are in cities. How can these be the richest areas for mammals?

Let’s have a look at the sampling effort (N: number of records per grid cell) and compare both layers

Mapping occurrence records 🦌

ggplot() + 
  geom_sf(data=mammalsCZ_grids, aes(fill=N)) +
  scale_fill_fermenter(palette ='YlGnBu', n.breaks=9, direction = 1) +
  geom_sf(data=CZ_borders, fill=NA) +
  coord_sf(crs=4326) +
  theme_bw()

Mapping occurrence records 🦌

What can you say about the hotspots of species richness we found?

We will see how to account for this issue later, while we wait for some geoprocessing steps :)

Mapping occurrence records 🦌

  1. Get occurrence records of mammals in the Czech Republic (POINTS)
  2. Create a grid layer of the Czech Republic (POLYGONS)
  3. Calculate sampling effort (N) and species richness per grid cell (SR)
  4. Visualise N and SR patterns in Czech Republic


We are done! No it’s your turn :)

GLOBAL Testudines of the World

Photo of Emys orbicularis observed in France by Martin Costechareire licensed CC-BY-NC, via iNaturalist

Mapping species ranges 🐢

Steps:

  1. Get testudines’ IUCN expert range maps (POLYGONS)
  2. Get a world map
  3. Create a global 1-degree grid (POLYGONS)
  4. Calculate the species richness per grid cell (SR)
  5. Visualise SR patterns in the world

Mapping species ranges 🐢

  1. Get testudines’ IUCN expert range maps

Read the data

testudines <- st_read('data/testudines.gpkg')
Reading layer `data_0' from data source 
  `/Users/flograttarola/Documents/GitHub/Spatial-Ecology-and-Macroecology/Practical_classes/Week2_gridding_and plotting/data/testudines.gpkg' 
  using driver `GPKG'
Simple feature collection with 169 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.999 ymin: -47.55406 xmax: 179.999 ymax: 61.50515
Geodetic CRS:  WGS 84


Polygons were downloaded from IUCN Spatial Data Download’s page

Mapping species ranges 🐢

  1. Get testudines’ IUCN expert range maps

What do the data look like?

testudines
Simple feature collection with 169 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.999 ymin: -47.55406 xmax: 179.999 ymax: 61.50515
Geodetic CRS:  WGS 84
First 10 features:
   ASSESSMENT ID_NO                 BINOMIAL PRESENCE ORIGIN SEASONAL
1      495907 10041      Heosemys annandalii        1      1        1
2      499158 10825    Indotestudo forstenii        1      1        1
3      499618 10950    Pangshura sylhetensis        1      1        1
4      505402 12124         Lissemys scutata        1      1        1
5      507698 12695      Malaclemys terrapin        1      1        1
6      507698 12695      Malaclemys terrapin        1      5        1
7      508210 12696   Malacochersus tornieri        1      1        1
8      508518 12775        Manouria impressa        1      1        1
9      511526 13038 Melanochelys tricarinata        1      1        1
10     511526 13038 Melanochelys tricarinata        6      1        1
        COMPILER YEAR                      CITATION
1  Anders Rhodin 2020 Chelonian Research Foundation
2  Anders Rhodin 2020 Chelonian Research Foundation
3  Anders Rhodin 2020 Chelonian Research Foundation
4  Anders Rhodin 2020 Chelonian Research Foundation
5  Anders Rhodin 2018 Chelonian Research Foundation
6  Anders Rhodin 2018 Chelonian Research Foundation
7  Anders Rhodin 2018 Chelonian Research Foundation
8  Anders Rhodin 2020 Chelonian Research Foundation
9  Anders Rhodin 2018 Chelonian Research Foundation
10 Anders Rhodin 2018 Chelonian Research Foundation
                                 LEGEND SUBSPECIES SUBPOP DIST_COMM ISLAND
1                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
2                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
3                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
4                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
5                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
6  Extant & Origin Uncertain (resident)       <NA>   <NA>      <NA>   <NA>
7                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
8                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
9                     Extant (resident)       <NA>   <NA>      <NA>   <NA>
10                   Presence Uncertain       <NA>   <NA>      <NA>   <NA>
   TAX_COMM                           geom
1      <NA> MULTIPOLYGON (((102.6212 5....
2      <NA> MULTIPOLYGON (((120.8806 1....
3      <NA> MULTIPOLYGON (((92.29536 22...
4      <NA> MULTIPOLYGON (((95.29987 15...
5      <NA> MULTIPOLYGON (((-81.45417 2...
6      <NA> MULTIPOLYGON (((-64.7296 32...
7      <NA> MULTIPOLYGON (((32.93152 -9...
8      <NA> MULTIPOLYGON (((101.684 5.0...
9      <NA> MULTIPOLYGON (((85.58097 25...
10     <NA> MULTIPOLYGON (((92 21.59583...

Mapping species ranges 🐢

  1. Get a world map

We will download a polygon of the world at medium scale and we will combine all countries into a unique polygon

world <- ne_countries(scale = 50, returnclass='sf')
world <- st_union(world) %>% st_make_valid() %>% st_cast() # fixes any problems


Let’s see how it looks

Mapping species ranges 🐢

  1. Get a world map
ggplot() + 
  geom_sf(data=world, fill='white')

Mapping species ranges 🐢

  1. Get a world map

Before working with this map, we need to project the layer: from lat/lon to equal area projection.


Earth is not flat :) Projections help us represent the two-dimensional curved surface of the globe into 2D space. There are many ways to do this. Find here a cool (Projection Wizard)[https://projectionwizard.org].

Mapping species ranges 🐢

  1. Get a world map

Equal-area maps preserve area measure, generally distorting shapes in order to do that

We will choose Equal Earth (EPSG:8857).

Mapping species ranges 🐢

  1. Get a world map

Let’s transform the data (both the world and the testudines’ layers)

world_ea <- st_transform(world, crs = 'EPSG:8857') %>% 
  st_make_valid() %>% st_cast()

testudines_ea <- st_transform(testudines, crs = 'EPSG:8857') %>% 
  st_make_valid() %>% st_cast()

Mapping species ranges 🐢

  1. Get a world map

Let’s double check that everything is alright

st_crs(world_ea, parameters = TRUE)$epsg
[1] 8857
st_crs(world_ea, parameters = TRUE)$ud_unit
1 [m]
st_crs(testudines_ea, parameters = TRUE)$epsg
[1] 8857
st_crs(testudines_ea, parameters = TRUE)$ud_unit
1 [m]

Mapping species ranges 🐢

  1. Create a global 1-degree grid (POLYGONS)

We will do this using the function st_make_grid()

st_make_grid(
  x,
  cellsize = c(diff(st_bbox(x)[c(1, 3)]), diff(st_bbox(x)[c(2, 4)]))/n,
  offset = st_bbox(x)[c("xmin", "ymin")],
  n = c(10, 10),
  crs = if (missing(x)) NA_crs_ else st_crs(x),
  what = "polygons",
  square = TRUE,
  flat_topped = FALSE
)

Mapping species ranges 🐢

  1. Create a global 1-degree grid (POLYGONS)

For the cellsize argument we will chose 1 degree, which is ~100km (=100,000m)

world_grid <- st_make_grid(world_ea, 
                           cellsize=100000,
                           square=FALSE) %>%  # this will make hexagons
  st_intersection(world_ea) %>% 
  st_sf(gridID=1:length(.))


We have the grid, we are ready to calculate metrics per grid cell. Let’s see how it looks

Mapping species ranges 🐢

  1. Create a global 1-degree grid (POLYGONS)

Mapping species ranges 🐢

  1. Calculate the species richness per grid cell (SR)

To do this, let’s keep only the testudines that are extant (resident) species

testudines_extant <- testudines_ea %>% 
  filter(LEGEND=='Extant (resident)') %>% 
  group_by(ID_NO) %>% 
  summarise(BINOMIAL=unique(BINOMIAL),
            PRESENCE=sum(PRESENCE)) %>% 
  ungroup() %>% st_cast()


This can take a while ☕️

Mapping species ranges 🐢

  1. Calculate the species richness per grid cell (SR)

Let’s plot 10 species as an example

Mapping species ranges 🐢

  1. Calculate the species richness per grid cell (SR)

Now let’s use st_join() to calculate the number of species per grid-cell (SR).

testudines_grid <- st_join(world_grid, testudines_extant) %>%
  group_by(gridID) %>%
  summarise(N=sum(!is.na(BINOMIAL)), 
            SR=n_distinct(BINOMIAL, na.rm = TRUE))


This will take a while ☕️🍪

Mapping species ranges 🐢

  1. Calculate the species richness per grid cell (SR)

Here are the results :)

testudines_grid
Simple feature collection with 20146 features and 3 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: -16920570 ymin: -8392928 xmax: 16920570 ymax: 8315130
Projected CRS: WGS 84 / Equal Earth Greenwich
# A tibble: 20,146 × 4
   gridID     N    SR                                                          .
    <int> <dbl> <int>                                             <GEOMETRY [m]>
 1      1     0     0                                 POINT (-16920565 -2061608)
 2      2     2     2 MULTIPOLYGON (((-16903934 -2109412, -16904928 -2110323, -…
 3      3     2     2 MULTIPOLYGON (((-16905520 -2108496, -16907015 -2103886, -…
 4      4     2     2 POLYGON ((-16782581 -2193350, -16785819 -2191050, -167869…
 5      5     3     3 POLYGON ((-16805308 -1830105, -16810425 -1830932, -168160…
 6      6     2     2 POLYGON ((-16783058 -2419039, -16784544 -2418515, -167857…
 7      7     2     2 MULTIPOLYGON (((-16779389 -2191507, -16775928 -2194312, -…
 8      8     2     2 MULTIPOLYGON (((-16749600 -2289778, -16749039 -2290664, -…
 9      9     2     2 POLYGON ((-16692191 -601861.4, -16691700 -601711.1, -1669…
10     10     2     2 POLYGON ((-16658328 -2432319, -16657130 -2432783, -166574…
# … with 20,136 more rows

Mapping species ranges 🐢

  1. Visualise SR patterns in the world

And now, we plot!

References

Any doubts?